/* Robustness ID */
clear all
cd "C:\Users\Public\Documents\ImmPanelRevis19\"
cd "ImmPanelRevis19\Donnees\select"
use final.dta, clear
gen reg=substr(ze,1,2)
gen nreg="11" if reg=="11" | reg=="00"
replace nreg="24" if reg=="24"
replace nreg="27" if reg=="26" | reg=="43"
replace nreg="28" if reg=="23" | reg=="25"
replace nreg="32" if reg=="31" | reg=="22"
replace nreg="44" if reg=="42" | reg=="21" | reg=="41"
replace nreg="52" if reg=="52" 
replace nreg="53" if reg=="53" 
replace nreg="75" if reg=="72" | reg=="74" | reg=="54"
replace nreg="76" if reg=="91" | reg=="73"
replace nreg="84" if reg=="83" | reg=="82"
replace nreg="93" if reg=="93" 
replace nreg="94" if reg=="94" 

xi , pre(REG_) i.nreg
xi , pre(CZ_) i.ze
save finalRobust.dta, replace

capture program drop autoresid
program define autoresid
version 11
syntax, qreg(name)
save tempX.dta, replace
predict resid , residuals
collapse (mean ) resid (count) nobs=resid, by(ze an)
egen indiv=group(ze)
egen time=group(an)
xtset indiv time
display "`qreg'"
gen dresid=D.resid
ivreg2 resid L.resid , cl(ze) noc
estimates store f`qreg'
ivreg2 resid L2.resid , cl(ze) noc
estimates store s`qreg'
use tempX.dta, replace
end

capture program drop robustID
program define robustID
version 11
syntax, occp(name) outc(name)

use finalrobust.dta, clear

/* create some outcome variables */
gen ldp=ln(dp/l.dp)
gen dftfy=(dp==360)-(L.dp==360)
gen oshift=(`occp'==0 & L.`occp'==1)

/* keep individuals in the occupation group */
keep if (btime==1 & `occp') | (btime==2 & L.`occp')

/* location in the decile of the initial group distribution */
gen lsnd=ln(sn/dp)
/* residual wages for each year */
quietly: reg lsnd F_* 
predict rlsnd, residuals
/* change in residual wages */
gen drlsnd=D.rlsnd
/* residual outcome */
quietly: reg `outc' F_*  if btime==2
predict r`outc' if btime==2, residuals

/* rank in the wage distribution */
drop F_*
egen qinit = fastxtile(rlsnd) if ftreat, by(ze an btime) nq(100)
replace qinit=qinit/100

/* interaction term and instruments */
gen lqinit=L.qinit
gen inter= dimm*L.qinit
gen linter=ldimm*L.qinit
gen inter75= dimm75 * L.qinit
gen linter75=ldimm75*L.qinit
gen intertm2= dimtm2 * L.qinit
gen lintertm2=ldimtm2 *L.qinit

keep if btime==2 & L.`occp' & !missing(dimm) & !missing(dimm75)
foreach year in 1982 1991 1999 2007 {
/* winsorize drlsnd */
quietly: sum drlsnd if an==`year', d
replace drlsnd=. if drlsnd>`r(p99)' & an==`year'
replace drlsnd=. if drlsnd<`r(p1)' & an==`year'
}
drop if missing(drlsnd)
drop if missing(basman)
drop if missing(bascom)

/* weight inverse of the size of the group in the CZ */
save temp.dta, replace
gen dummy=1
collapse (sum) nbze = dummy , by(ze an)
gen iw2=1/ nbze
drop nbze

save iw2.dta, replace
use temp.dta, replace
capture drop _merge
joinby ze an using iw2.dta , unm(m)
drop _merge

save temp1.dta, replace

use temp1.dta, clear
/* instrument baseline 75 */
ivreg2 r`outc' (dimm inter = dimm75 inter75) lqinit  Y_*  REG_* [aweight = iw2], /*partial(Y_* ) */ cl(ze)
estimates store `occp'75
autoresid , qreg(`outc'x75)

/* instrument baseline tm2 */
ivreg2 r`outc' (dimm inter = dimtm2 intertm2) lqinit  Y_*  REG_* [aweight = iw2], /*partial(Y_* ) */  cl(ze)
estimates store `occp'tm2
autoresid , qreg(`outc'tm2)

/* include lag */
ivreg2 r`outc' (dimm ldimm inter linter = dimtm2 intertm2 ldimtm2 lintertm2) lqinit  Y_*  [aweight = iw2], partial(Y_* ) cl(ze)
estimates store `occp'lag
autoresid , qreg(`outc'lag)

/* CZ fixed effect */
ivreg2 r`outc' (dimm  inter  = dimm75 inter75) lqinit  Y_* CZ_* [aweight = iw2], partial(Y_* CZ_*) cl(ze)
estimates store `occp'czfe
autoresid , qreg(`outc'czefe)

end

capture rm "C:\Users\Public\Documents\Tab21\Table13.rtf"
/* OUTFLOWS */
robustID, occp(bc) outc(lshift)

/* IVSELECT */
esttab bc75 bctm2 bclag bcczfe ///
using "C:\Users\Public\Documents\Tab21\Table13.rtf" /// 
, append title("lshift")  b(%9.3f) cells(b(star fmt(3)) se(fmt(3) par) ) ///
 stats(N widstat) star(* 0.10 ** 0.05 *** 0.01)
 
/* autocorrelation tests */
estout flshiftx75 flshifttm2 flshiftlag  flshiftczefe  , cells(b(star fmt(%9.3f)) ///
se(par(`"="("'`")""')) t p ) stats(N r2 widstat) starlevels(* 0.10 ** 0.05 *** 0.01)
estout slshiftx75 slshifttm2 slshiftlag slshiftczefe  , cells(b(star fmt(%9.3f)) ///
se(par(`"="("'`")""')) t p ) stats(N r2 widstat) starlevels(* 0.10 ** 0.05 *** 0.01)

/*************************************/
/* Inflows into the occupation group */
/*************************************/
capture program drop inflowsrobid
program define inflowsrobid
version 11
syntax, occp(name)

use unbalanced.dta, clear
/*gen cad2=(cs2h=="33" | cs2h=="34" | cs2h=="37" | cs2h=="38")*/

gen orig=(ptreat==1 & `occp'==1 )
/* gen inflow=(dtreat==1 & bc==1 & (L.bc!=1 | L.ize!=ize) & !missing(L.ize)) */
gen inflow=(dtreat==1 & `occp'==1 & L.ize!=ize )


gen reg=substr(ze,1,2)
gen nreg="11" if reg=="11" | reg=="00"
replace nreg="24" if reg=="24"
replace nreg="27" if reg=="26" | reg=="43"
replace nreg="28" if reg=="23" | reg=="25"
replace nreg="32" if reg=="31" | reg=="22"
replace nreg="44" if reg=="42" | reg=="21" | reg=="41"
replace nreg="52" if reg=="52" 
replace nreg="53" if reg=="53" 
replace nreg="75" if reg=="72" | reg=="74" | reg=="54"
replace nreg="76" if reg=="91" | reg=="73"
replace nreg="84" if reg=="83" | reg=="82"
replace nreg="93" if reg=="93" 
replace nreg="94" if reg=="94" 

gen sond=12
replace sond=24 if an==2007
destring nreg, replace
destring ze, replace

save temp1.dta, replace
collapse (sum) norig=orig ninflow=inflow [fweight = sond], by(ze an)
save p1.dta, replace
use temp1.dta, clear


collapse dimm dimm75 ldimm ldimm75 dimtm2 ldimtm2 nreg Y_*, by(ze an)
joinby ze an using p1.dta, unm(m)
gen poid=sqrt(norig)

drop _merge

egen ize=group(ze)
egen time=group(an)
xtset ize time
gen inflr=(D.ninflow)/norig

xi , pre(REG_) i.nreg
xi , pre(CZ_) i.ze

keep if !missing(inflr)
drop time
/* test for serial correlation of residuals */
egen indiv=group(ze)
egen time=group(an)
xtset indiv time

/* IV */
ivreg2 inflr (dimm  = dimm75 ) Y_*  if !missing(ldimm75), partial(Y_* ) cl(ze)
estimates store `occp'75
abar , l(1)
abar , l(2)

ivreg2 inflr (dimm  = dimtm2 ) Y_*  if !missing(ldimm75), partial(Y_* ) cl(ze)
estimates store `occp'tm2
abar , l(1)
abar , l(2)

ivreg2 inflr (dimm  ldimm = dimm75 ldimm75 ) Y_* , partial(Y_* ) cl(ze)
estimates store `occp'lag
abar , l(1)
abar , l(2)

ivreg2 inflr (dimm   = dimm75  ) Y_* REG_*, partial(Y_* REG_*) cl(ze)
estimates store `occp'rfe
abar , l(1)
abar , l(2)

ivreg2 inflr (dimm   = dimm75  ) Y_* CZ_*, partial(Y_* CZ_*) cl(ze)
estimates store `occp'czfe
abar , l(1)
abar , l(2)


end

inflowsrobid, occp(bc) 

/* IVSELECT */
esttab bc75 bctm2 bclag bcczfe ///
using "C:\Users\Public\Documents\Tab21\Table13.rtf" /// 
, append title("INFLOWS")  b(%9.3f) cells(b(star fmt(3)) se(fmt(3) par) ) ///
 stats(N widstat) star(* 0.10 ** 0.05 *** 0.01)
 
/*********************************/
/* WAGES WAGES WAGES WAGES WAGES */
/*********************************/
capture program drop wagesrobid
program define wagesrobid
version 11
syntax, occp(name) 

/********************/
/* Balanced Sample  */
/********************/
use finalRobust.dta, clear

/* keep individuals in the occupation group */
keep if (btime==1 & `occp') | (btime==2 & L.`occp')

/* location in the decile of the initial group distribution */
gen lsnd=ln(sn/dp)
/* residual wages for each year */
quietly: reg lsnd F_* 
predict rlsnd, residuals
/* change in residual wages */
gen drlsnd=D.rlsnd

gen ldp=ln(dp/l.dp)
gen oshift=(`occp'==0 & L.`occp'==1)
keep if btime==2 & L.`occp' & !missing(dimm) & !missing(dimm75)
foreach year in 1982 1991 1999 2007 {
/* winsorize drlsnd */
quietly: sum drlsnd if an==`year', d
replace drlsnd=. if drlsnd>`r(p99)' & an==`year'
replace drlsnd=. if drlsnd<`r(p1)' & an==`year'
}
drop if missing(drlsnd)
drop if missing(basman)
drop if missing(bascom)

/* weight inverse of the size of the group in the CZ */
save temp.dta, replace
gen dummy=1
collapse (sum) nbze = dummy , by(ze an)
gen iw2=1/ nbze
gen poid=sqrt(nbze)

save iw2.dta, replace
use temp.dta, replace
capture drop _merge
joinby ze an using iw2.dta , unm(m)
drop _merge

save temp1.dta, replace
collapse drlsnd dimm dimm75 dimtm2 ldimm ldimm75 ldimtm2 Y_* REG_* CZ_* poid, by(ze an)

/* test for serial correlation of residuals */
egen indiv=group(ze)
egen time=group(an)
xtset indiv time

/* IV */
ivreg2 drlsnd (dimm  = dimm75 ) Y_* , partial(Y_* ) cl(ze)
estimates store `occp'75
abar , l(1)
abar , l(2)

/* IV */
ivreg2 drlsnd (dimm  = dimtm2 ) Y_* , partial(Y_* ) cl(ze)
estimates store `occp'tm2
abar , l(1)
abar , l(2)

/* IV */
ivreg2 drlsnd (dimm ldimm = dimm75 ldimm75 ) Y_* , partial(Y_* ) cl(ze)
estimates store `occp'lag
abar , l(1)
abar , l(2)

/* CZ FE */
ivreg2 drlsnd (dimm  = dimm75 ) Y_* CZ_* , partial(Y_* CZ_*) cl(ze)
estimates store `occp'czfe
abar , l(1)
abar , l(2)

end

/* WAGES */
wagesrobid, occp(bc) 
/* IVSELECT */
esttab bc75 bctm2 bclag bcczfe ///
using "C:\Users\Public\Documents\Tab21\Table13.rtf" /// 
, append title("wages")  b(%9.3f) cells(b(star fmt(3)) se(fmt(3) par) ) ///
 stats(N widstat) star(* 0.10 ** 0.05 *** 0.01)

/* other groups pour stocker !! */
 